Hull 


AD“A257  941 

■mtiiiii 


<EKXL.T3;  vers  1;  0V24/9j> 


EKXL:  a  Dynamic  Poisson-Solver 


Katharine  E.  King  and  Robert  W.Bussard 


SDTIC 

A  U 


EMC2H)791-03 


WROVEDFORr.'a.'CRaE 

£>-<W3'r;ovu-M,v:,i;; 


^EASE 
TED  ■ 


this  document  has  been  approved 
for  public  release  ond  sole;  its 
distention  is  unlimited. 


CiE/lHi-rt 

OCT  2  2  1992  "  4 


M.'  r., . 

......  .  CF  ' ' 


the  . l 


||jg9994 


iMC 


rACIBAl  ACCBSacS  08 


=srK.%“s-»s«» 


^'7/ 


1. 


fflSTORY  AND  BACKGROUND  OF  THE  EKXL  CODE 


The  EKXL  code,  and  each  of  its  predecessors,  is  a  1— D  radially-dependent  Poisson-solver. 
It  gives  static  solutions  for  potential  and  density  distributions,  given  edge  density 
specifications  for  these  parameters.  It  is  a  "point"  model  in  that  spatial  variation  of 
externally-imposed  fields  is  possible,  but  spatial  dependence  of  inter— particle  interaction 
phenomena  can  not  be  accomodated  as  an  integral  part  of  code  operations. 

The  original  code,  IACCEL2/EACCEL2,  was  developed  by  Ensleyi-  ^  and  provided  closed 
analytic  solutions  in  integral  form  using  Struve  and  Bessel  functions.  This  code  employed 
Gaussian  distributions  for  transverse  energy  and  momentum.  Solutions  were  limited, 
however,  to  electron  currents  of  1—10  A  and  core  densities  of  1x10^—1x10^®  1/cm^,  as 
indicated  in  Figure  1. 


The  first  modification  to  Ensley’s  work  was  ACCEL,  developed  by  Mission  Research 
Corporation^  (MRC).  This  version  of  the  code  utilized  a  more  efficient  algorithm  for 
converging  a  solution  to  Poisson’s  equation,  resulting  in  significant  run-time  savings.  New 
physics  elements  were  also  added  in  order  to  indude  tangential  current  densities  where 
previously  only  radial  current  was  used.  As  shown  in  Figure  1,  with  the  ACCEL 
improvements  solutions  were  achieved  with  electron  currents  of  100-500  A  and  core 


densities  of  1x10^^— 1x10*^  1/cm*. 

What  followed  was  a  new  code,  XL,  also  developed  by  MRC<,  that  was  utilized  to  its 
capadty  by  EMC^  to  study  the  Polywdl*“/HEPS  system*.  This  code  employed  a 
differential  numerical  analysis  procedure  and  matrix  inversion  solution  method.  The 
distribution  functions  for  transverse  energy  and  momentum  were  changed  to  square 
functions  to  simplify  the  calculation.  Features  added  at  this  time  were  adjustable  precision 
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and  an  adaptive  grid.  This  code  achieved  much  higher  electron  recirculating  currents, 

1x10  A,  and  core  densities,  1x10  -1x10  1/cm  .  However,  the  core  density  barrier  at 
1x10^^  1/cm^  could  not  be  surpassed,  making  any  studies  of  Polywell^”'  reactor  densities 
impossible.  Work  was  then  underway  at  MRC  on  the  next  revision,  KXL. 

n.  THE  ADIABATIC  KXL  CODE 

O 

EMC  began  working  with  the  adiabatic  code  KXL  in  July,  1990.  KXL  uses  the  same 
mathematical  engine  developed  in  the  XL  code  and  a  series  of  intermediate  "steady-state" 
solutions  over  time  to  reach  the  final  steady-state  condition  for  the  given  input 
parameters.  Selected  parameters,  such  as  input  current  and  the  length  of  the  time  step, 
can  be  adjusted  as  the  calculation  proceeds.  As  indicated  in  Figure  1,  this  method  and 
some  of  the  improvements  discussed  below  allow  the  code  to  reach  much  higher  core 
densities  than  were  possible  with  XL. 

The  original  MRC  version  of  KXL  used  an  explicit  input  of  confinement  time.  This  was 
replaced  with  input  of  the  recirculation,  G.  j,  and  the  confinement  time  is  calculated  in  the 
code  by  G.  .t^  .  This  allowed  for  more  consistency  since  the  transit  time  is  also  used 

elsewhere  in  the  code. 

The  "trapping"  formula,  that  allows  feedback  on  the  current  density  was  first 
changed  by  David  Smithe  of^piC  and  then  later  by  EMC2  (see  bdow).  The  original  form 
was  a  generic  adjustment  that  reduced  the  input  current  by  an  input  reduction  factor. 

This  was  changed  to  an  algorithm  that  adjusts  the  current  based  on  the  core  potential  to 
prevent  the  anode  from  maximizing,  due  to  too  many  ions  in  the  core,  causing  the 
calculation  to  be  unsuccessful.  A  target  anode  potential  and  maximum  beam  current  are 
input  and  the  code  adjusts  the  beam  current  up  to  that  maximum  in  an  attempt  to  obtain 
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the  target  anode  height.  The  input  variable  is  and  it  represents  the  target  well  depth. 
Vg  can  be  different  for  each  species,  allowing  separate  control  over  each  beam.  For 
example,  a  for  the  electron  beam  greater  than  the  maximum  well  depth  would  set 
equal  to  one  and  allow  the  electron  beam  to  remain  on  at  maximum  current  for  the  entire 
run  while  the  code  adjusted  the  ion  current  within  the  indicated  range.  The  anode  control 
formula  is  as  follows: 


F*  = 

trap 


where  a  = 


l+e‘ 


^.03  <p 


(1) 


max 


(p^  =  Current  Core  Potential,  ip  =  Maximum  Core  Potential,  and  the  *  is 
^<1  '  ^max  * 


determined  by  the  charge  on  the  beam. 


If  the  maximum  available  ion  current  is  too  high,  however,  the  F^^^^  formula  cannot  turn 
the  current  down  enough.  In  effect,  there  is  gun  "leakage"  and  the  ions  build  up  (due  to 
their  recirculation)  and  overwhelm  the  core.  When  using  a  large  ion  current,  it  still  needs 
to  be  within  a  reasonable  range  of  the  final  gun  current  for  the  code  to  run  successfully. 


The  auxiliary  potential  was  introduced  in  KXL  to  model  mirror  reflection  of  the  electrons 
in  the  cusp  region,  /iVB  ==  jmv^.  The  input  to  the  code  is  the  magnitude  of  therefore, 
it  is  not  evident  what  is  the  magnitude  of  B  since  /x  is  defined  as  a  function  of  both  B  and 
v^.  Much  analysis  went  into  the  auxiliary  potential  in  an  attempt  to  determine  the 
magnitude  of  B  to  use  in  the  ^culation  based  on  either  the  wiffle  ball  or  mirror 
reflection  modds,  svntching  between  the  two  where  appropriate.  Since  the  magnitude  of  B 
is  indeterminate  from  the  auxiliary  potential,  this  linkage  was  not  possible. 

In  previous  work  with  the  XL  code  a  gridding  proUem  had  been  identified.  From  extensive 
code  runs,  many  of  which  did  not  succeed,  it  was  found  that  as  the  system  density 
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increased  the  adaptive  grid  was  not  correctly  repositioning  the  grid  points  to  allow  proper 
resolution  of  areas  of  high  particle  density.  A  radial  weighting  factor  ("R"  factor)  was  then 
used  to  ensure  an  adequate  number  of  grid  points  near  the  outer  radius,  the  region  where 
the  density  changes  rapidly.  This  weighting  of  the  outer  regions  did  not  allow  sufficient 
grid  points  in  the  core  region  to  resolve  high  density  cases.  After  EMC^  and  MRC  had 
worked  out  a  solution  for  the  XL  code,  the  "R”  factor  was  reintroduced  when  the  KXL 
code  was  written.  Once  this  factor  was  reinoved  from  KXL  the  cases  began  to  run  more 
consistently  and  the  core  resolution  was  greatly  improved.  Figure  2  a  shows  a  KXL  run 
with  the  "R"  factor  and  Figure  2  b  shows  the  same  case  after  the  "R"  factor  was  removed. 
Note  the  improved  resolution  of  the  core  region  in  Figure  2  b 


MRC  compiled  the  KXL  code  on  the  Cray  computer,  and  previous  PC  cases  were  run  to 
verify  that  both  codes  were  consistent.  Once  the  running  parameters  for  ion  current  were 
worked  out,  Cray  cases  were  run  to  reach  the  high  densities  that  appeared  unattainable  on 
the  PC,  in  an  attempt  to  fill  in  high  density  points  in  a  study  of  core  density  vs.  electron 
current.  Examining  this  data,  however,  it  became  apparent  that  the  core  densities  were 
much  lower  than  expected  from  earlier  calculations.  Numerous  variations  of  input 
variables  were  tried  in  an  attempt  to  increase  the  densities.  In  the  mean  time,  an 
explanation  was  sought  as  to  why  the  densities  w^e  lower  than  predicted,  and  attempts 
were  made  to  model  this  behavior.  Also,  the  slope  of  core  density  vs.  electron  current  was 
expected  to  be  linear  for  constant  electron  confinement.  As  Figure  3  indicates,  it  was  not 
linear  until  the  electron  curre%  reaches  1000  A.  This  problem  was  discussed  at  some 
length  by  EMC^,  MRC  and  Krall  Associates. 

The  explanation  for  the  change  in  slope  at  1000  A  was  found,  but  without  illuminating  the 
too-low  density  problem.  The  dectron  density  is  increasing  linearly  the  entire  time,  as 
expected.  It  is  only  the  ion  density  that  changes,  due  to  the  anode  height.  The  fn  required 
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A  A 

to  form  an  anode  is  appioximatdy  &i  =  (n.  -  n^)  :=  3.32x10  The  electron 

currents  below  the  linear  slope  corresponded  to  densities  where  &i  was  a  significant  portion 
of  the  core  density.  As  core  density  moves  above  this  density,  6n  becomes  less  significant, 
the  effect  disappears,  and  the  slope  finally  becomes  linear.  By  plotting  this  data  for 
varying  anode  heights  it  is  apparent  that  the  curve  approaches  the  linear  condition  as  the 
anode  height  approaches  zero  (Figure  4). 


In  an  effort  to  understand  the  too-low  core  densities,  attention  turned  to  the  transit  times 
The  code  was  modified  to  print  the  transit  time  at  each  step  and  this  information  was 
plotted,  revealing  that  the  code-calculated  transit  times  were  much  too  small  to 
correspond  to  reality  (Figure  5).  Several  changes  were  suggested  and  tried  to  improve  the 
transit  times,  however  it  was  then  found  that  altering  the  transit  time  calculation  had  no 
effect  on  the  final  solution.  This  led  to  a  suspicion  that  the  dn/dt  equation  that  was  used 
to  link  the  intermediate  steady-state  ("snapshot”)  solutions  might  be  incorrect.  In  the 
KXL  code  this  equation  (for  electrons)  had  the  form: 


dt 


(2a) 


which  gives  the  final  steady-state  relation  =  n^.  The  basic  premise  of  this 
time-dependent  coupling  equation  (2a)  is,  in  fact,  false.  The  correct  formulation  is  one 
that  coimects  time-depende^  of  total  electron  populations  in  the  system  between 
successive  static  solutions  from  the  KXL  code.  This  equation  is: 


dn 

dt 


2  « 
4irR  n.  .V.  .  n 
_ »°J _ 


4/3  tR*  t^ 


(2b) 
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However,  since  the  codes  are  all  written  to  solve  only  one-dimensional  boundary-value 
problems,  it  was  necessary  to  recast  this  equation  (2b}  in  terms  of  edge  (boundary) 
densities.  Thus  major  modifications  to  this  portion  of  the  code  were  required.  And,  other 
physics  phenomena  of  importance  to  system  operation  also  needed  to  be  added  to  the  code 
to  improve  the  realism  of  its  output  results. 

Taking  all  of  these  factors  together,  at  this  point  it  was  apparent  that  the  KXL  code 
needed  considerable  modification  in  order  to  be  useful  in  analyses  of  Polywell^"'  machines 
in  this  (HEPS)  program.  Leaving  the  original  generic  code  intact,  EMC  then  began 
development  of  the  EKXL  code,  using  the  mathematical  engine  of  KXL  as  developed  by 
MRC  (and  EMC^)  but  now  including  all  of  the  appropriate  ion/electron/plasma  physics 
analytic  forms  appropriate  to  and  required  for  correct  numerical  results  from  EKXL  code 
computions. 

ra.  THE  EKXL  CODE 


Change  in  Edge  Density  with  Respect  to  Time 


The  first  step  in  customizing  EKXL  to  the  PolyweU**”  system  was  to  derive  a  new,  correct 
dn/dt  equation.  This  equation  is  the  key  to  the  linkage  of  intermediate  steady-state 
solutions  over  time,  and  numerical  tests  have  shown  solutions  to  be  very  sensitive  to  errors 
in  this  equation.  The  change  the  total  number  of  particles  in  the  system  is  simply  the 
input  less  the  losses  over  a  given  period  of  time.  Therefore,  the  change  in  the  average 
number  of  particles  with  respect  to  time  is  given  by: 


dn  4jrR^n.  .v.  .  n  3n.  .  n 

_ _ _ inj  inj _ _  iBj _ 

dt”  4/3  tR’  R/Vj.j  t. 


(3) 
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where  is  the  confinement  time,  G.  The  code  is  only  a  differencing  solution  of  a 

one-dimensional  boundary-value  problem.  Thus  dn/dt  can  only  be  tied  to  the  change  in 
edge  (boundary)  density,  as  opposed  to  average  density,  so  a  parameter,  =  n/n^,  was 
introduced  to  relate  edge  density  to  average  density.  Since  this  ratio  is  not  fixed,  is 
calculated  in  the  code  each  time  the  dn/dt  equation  is  solved  (see  below).  A  second 
parameter,  the  ratio  of  free-stream  half-transit  time  with  no  potential  well  to  the  actual 
transit  time. 


R/v.  . 

p  - 

“  t ,  ’ 

trans 


(4) 


was  introduced  to  allow  the  transit  time  to  be  factored  out  of  the  dn/dt  equation.  The 
final  form  of  the  dn/dt  equation  is  then: 


dn 

e 

dt 


F  F 

n  t 


t 


trans 


Calculation  of  F 

n 


(5) 


The  average  number  of  particles  between  each  pair  of  grid  points  is  used  to  calculate  the 
number  of  particles  in  that  shell.  This  is  summed  over  the  entire  grid  and  divided  by  the 
total  volume  to  arrive  at  n.  A  modification  was  necessary,  however,  for  the  electron 
calculation  in  the  outer  shdl.*^  high  density  cases,  the  dectron  density  diange  across  the 

«  j 

last  grid  space  can  be  very  large  (ca.  10  — 10  x)  and  an  algebraic  average  is  not  accurate. 
An  exponential  mean  form  was  derived  based  on  the  assumption  that  the  logrithmic 
gradient  of  the  edge  electron  density  varies  inversely  with  the  Debye  length,  at  the 
outer  radius,  R,  d/dr  (In  n^)  =  with  a  parameter,  fi,  (typically  0  similar  to 
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0.7),  that  characterizes  the  drop  in  E  between  the  outer  radius  and  the  first  grid  point, 
<rj>.  When  the  electron  density  begins  to  increase  rapidly  in  this  region,  EKXL  will 
switch  from  an  algebraic  average  to  the  following: 


(I-<r,>)  R 


Electron  and  Ion  Recirculation 


Although  analytic  forms  have  been  derived  for  ion  loss  due  to  core  upscatteringB*  ^  and 
from  fusion  reactions*,  it  has  been  found  convenient  to  choose  a  constant  ion  recirculation, 
G.,  as  an  independent  input  to  the  code.  Variation  in  G.  affects  only  the  ion  injection 
current,  I.,  as  the  circulating  ion  current,  I  =  I.G.,  is  fixed  by  the  parameters  of  a  given 
case.  The  code  has  the  ability  to  adjust  I.  as  needed  (see  below). 

The  electron  recirculation,  Gj,  may  be  input  at  a  constant  value  or  it  may  be  calculated  in 
the  code  using  the  physics  of  electron  confinement  under  mirror  reflection  (MR)  and  wiffle 
ball  (WB)  modes.  This  model  is  discussed  in  detail  in  other  EMC2  technical  notes* <  and 
its  pertinent  algorithms  and  defining  equations  are  given  in  Figure  6.  As  discussed  earlier, 
the  magnitude  of  B  cannot  be  tied  to  the  auxiliary  potential  and  is,  therefore,  a  separate 
input  parameter  for  the  mirror/wiffle  model. 

Ion  Losses  Due  to  Fusion  Reactions 


Fusion  reactions  of  ions  in  the  core  and  surrounding  region  will  consume  ions.  This  source 
of  loss  is  important  only  at  conditions  of  high  core  density  and  large  ion  energy.  The  total 
fusion  rate  is  calculated  in  the  code  as  part  of  the  fusion  and  gain  calculations  (see  below). 
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This  allows  for  easy  inclusion  of  the  ion  losses  due  to  fusion  in  the  ion  dn/dt  equation.  If 
the  total  fusion  rate  is  fusions/sec,  then  the  number  of  ions  lost  to  fusion  is 
ions/sec.  This  loss  term  would  go  in  to  the  original  dN /dt  equation  and  the 
resulting  dn  /dt  equation  would  be: 


dn 
_ e 

dt 


f3n.  . 

‘“1 

’  n 

DD,  1 

-S  + 

f  u  • 

F  F, 

T 

G. 

F  Fr^T^ 

D  t 

1 

n  DD 

trans 


(7) 


where  =  4/3  /t^^^^^  is  introduced  to  allow  the  transit  time  to  be  factored  out. 

Ion- Losses  Due  to  Core  Upscattering 


Ion  losses  due  to  core  collisional  upscattering  can  be  accounted  for  by  adding  a  loss  term  to 
the  dn/dt  equation  for  ions.  Based  on  the  work  of  Rosenberg  and  Krall>,  as  modified  by 
Bussard^,  the  upscattering  time  for  ions,  t^^  is  calculated  each  time  the  dn/dt  equation  is 
solved,  from 

1  1 
t^p  1.892x10^^  f^(l+f)^  R 

A  factoring  parameter,  introduced  into  the  dn/dt  equation  so  that  the 

ion  losses  due  to  upscattering^e  given  by  n^/(F^t^^^^^).  This  form  is  easily  incorporated 
into  the  ion  dn/dt  equation  which  becomes 


n(r) 


E(r) 


3/2 


dr 


(8) 


dn 

e  __ 

'  3n.  , 

n 

ri  ‘  1 

DD,  1 

f  u  ■ 

1 

dt 

F  F,  * 

- 1 

G. 

F 

F  F„„ 

t. 

n  t 

W 

1 

u 

n  DD  ^ 

trani 

(9) 
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The  fractional  velocity  spread,  f,  is  a  specified  input  parameter.  If  f  =  0,  no  ion 
npscattering  is  calculated.  This  parameter  can  be  estimated  by  comparison  of  scattering 
collision  rates  in  the  core,  edge,  and  mantle  regions of  the  system,  comparing 
transverse  momentum  vs.  radial  momentum  exchange  processes.  Such  estimates  are  not 
yet  complete;  when  available  the  analytic  forms  that  describe  them  can  be  added  to  the 
dn/dt  equation,  as  above. 

Average  Transit  Time  Calculations 

The  data  in  the  code  is  quite  sufficient  to  calculate  the  ion  transit  time  directly.  Using  the 
potential  at  each  grid  point  and  the  input  parameters  £_,  dE  and  average 

transit  time  across  each  grid  point  can  be  calculated  and  summed  across  the  entire  grid  to 
arrive  at  the  average  ion  transit  time  through  the  system. 

The  calculation  of  transit  time  for  electrons,  however,  presents  an  interesting  problem. 

The  direct  method  used  for  ion  transit  time  does  not  lend  itself  to  an  accurate  calculation 

of  the  average  electron  transit  time  as  the  individual  times  vary  greatly  with  the  energy 

spread  of  the  dectrons.  Electrons  with  a  maximum  dE  will  be  turned  early  and  have 

perp 

short  transit  times  while  the  electrons  with  maximum  energy,  E.  .+dE,  and  no  dE  _ 

ujj  perp 

will  pass  through  the  bottom  of  the  well  and  reach  the  core  yielding  much  longer  transit 
times. 

Several  attempts  were  made  to  calculate  and  average  the  electron  transit  times  within  the 
code,  consuming  significant  computational  time.  Due  to  the  dependence  of  on  transit 
time,  however,  a  small  deviation  in  average  transit  time  could  continue  to  feed  back  into  a 
large  error  in  dn/dt.  In  short,  the  detailed  computational  capability  of  the  code  would  lead 
to  amplification  of  numerical  errors  in  original  data.  In  this  circumstance,  the  internal 
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data  calculation  was  abandoned  in  favor  of  an  analytic  form  based  upon  system 
parameters,  which  yields  a  more  accurate  and  consistent  averaged  result  and  utilizes  less 
computational  time.  As  derived  and  discussed  in  another  EMC2  note^^  (see  for  definitions 
of  terms)  the  electron  transit  time  is: 


_  2R 

* 

<ro>+ 

r<r  >^-<r  >^1 

w  0 

2<I,> 

V. 

1  n 

*  • 

1/2 

<r  > 

W 


1/2 


where  a  =  ip  Itp  ri'  =  1-a  (l-w  ),  n  =  E  /E  , 

q  ^o'  ^min’  '  q'  V’  'o  ao'  o’ 

r  2  2  1  l/(ro+2) 

<r  >  =  <r.>  +  —  <r.>  ri'  and  <r„>  is  the  core  radius, 

w  0  Im  0  '  I  0 


This  is  a  calculation  of  the  maximum  transit  time.  But  the  electrons  occupy  a  spread  in 
total  and  transverse  energy.  Across  this  energy  distribution  the  individual  transit  times 
can  vary  markedly.  Analysis  of  this  variation  suggests  an  exponential  distribution  of 
transit  times  across  the  band.  Integrated  averaging  over  this  variation  yields  a  correction 
factor  to  give  average  transit  times  of  the  average  electron  in  terms  of  the  total  maximum 
transit  time  of  equation  (10).  This  factor  is  found  to  be: 


In 


r  = 

trans 


t  rsns 


2K7?: 


in 


t  rant 


5K7VT 


-  1 


tot 


(11) 


II 


Fndoii  Rates  and  Power  Gain 


The  KXL  code  contained  a  subroutine  to  calculate  the  approximate  neutron  rate  (1/2  the 
fusion  rate)  for  DD  fusion.  This  subroutine  was  modified  in  EKXL  to  include  fusion  rate 
calculations  for  DT,  D^He  and  p^^B  based  upon  the  deuteron  density  (as  the  code  does  not 
currently  accommodate  different  electron  and  ion  injection  radii,  necessary  for 
center-of— momentum,  center-of-geometry  placement  of  these  ions).  Calculation  of 
system  gross  gain  for  all  four  fuel  combinations  was  also  added  to  the  EKXL  code,  where 
gross  gain  is  defined  as 


GGAIN  = 


(r)a(E)v.(r)47rr^dr 


P  +  P. 

e.  .  1.  . 

inj  inj 


(12) 


The  densities  used  for  each  fuel  combination  were: 

DT  For  DT  fusion  rates,  the  existing  deuteron  density,  D,  is  divided  in  half  and 
treated  as  50:50  D:T 

D^He  For  a  50:50  D:^He  distribution  the  total  D  +  ^He  density  will  be  2/3  the 
existing  D  density  for  the  same  potential  distribution. 

p^^B  For  50:50  p:^^^he  total  density  of  p  +  ^^B  is  1/3  the  existing  D  density. 


Bremmstrahlung  Losses 


Bre^mstrahlung  will  be  an  important  energy  loss  mechanism  from  the  system  at  high 
density  conditions.  This  was  estimated  from  standard  formulae's  for  bremmstrahlung 
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radiation  density,  integrated  over  the  system  volume  to  give  total  biemmstrahlung  output, 
ignoring  self-absorption  by  the  plasma  and  any  internal  reflection  from  the  boundary 
walls.  This  is  discussed  in  a  separate  EMC2  technical  note,  forthcoming.  The  resulting 
expression  is  used  in  the  EKXL  code  to  calculate  bremmstrahlung  gross  output  as 

Pbrenm.  = 

This  is  added  to  the  injection  power  to  determine  net  power  input  and  system  net  gain. 


R 


ni(r)n^(r)Z^VE^(r)47rT^dr 


(13) 


NGAIN  = 


fus 


p.  .  +p. 

mj  bremm 


Ion  Current  Adjustment 


(14) 


The  problem  with  the  anode  control,  discussed  previously  -  allowing  too  many  ions  to 
build  up  in  the  system  and  thus  requiring  the  maximum  ion  current  to  be  within  range  of 
the  final  ion  current  —  was  addressed  in  EKXL  with  a  different  form  of  anode  control  that 
calculates  the  ion  current  density  based  on  the  core  densities,  the  electron  current  density 
and  G.^.  The  difference  in  core  densities  is  given  by^* 


n.  —  n  =  3.32  x  10 

1  e 


>  (vO 


(15) 


where  is  the  core  potential  and  is  the  potential  at  the  bottom  of  the  well  (both  in 
eV),  r^  is  the  core  radius  in  meters  and  density  is  in  1/m  .  The  difference  in  densities  at 
the  target  anode  potential,  ip  ,  is  then 

ft 
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(16) 


a  (ip  -Ip  .  ) 

-  n'  =  3.32  x  10®  * 

I  e  _2 


The  change  in  density  from  n.  to  n  j  is  then 


o{ip-ip  ) 

An.  =  ni -n.  =  3.32x10® (An  +n)-n. 

iii  2  'ec'i 


or 


a  (iP  -IP  .  ) 

An.  =  3.32  X  10®  —  (n.  -n  )  +  An 

1  _2  '  1  e'  e 


(17) 


Substituting  equation  (15)  into  equation  (17)  gives 


a(ip  -ip) 

An.  «  3.32  X  10®  ■  *  ”  +  An 
1  ••  ® 


(18) 


and 


3.32  X  10®  G. 

*^i  = - 2 —  ( 

1  •o 


G. 


(19) 


where  J.  ^  are  the  circulating  current  densities. 


Automatic  Operating  Mode 


In  the  XL  code  a  case  achieved  a  steady-state  solution  or  none  at  all.  The  KXL  code  has 
the  ability  to  stop  at  any  "snapshot"  point  by  specifying  the  time  duration  in  the  input  file. 
This  allows  a  succeeding  run  to  be  started  at  the  end  of  a  previous  calculation,  using  a 
special  data  file  created  when  that  calculation  was  completed,  and  to  continue  for  the  time 
duration  specified  in  the  new  input  file.  This  allows  for  the  adjustment  of  time  steps  and 
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ion  current  necessary  to  achieve  the  final  steady-state  solutions  previously  unattainable 

with  the  XL  code.  The  drawback  is  that  the  manual  adjustment  of  duration  time  and  ion 

currents  require  continuing  operator  input.  In  the  EKXL  code  an  "automatic"  mode  was 

introduced  that  eliminates  the  need  for  most  manual  restarts.  In  automatic  mode,  the  code 

will  run  to  the  end  of  the  specified  time  duration  or  a  steady-state  solution,  whichever 

comes  first,  using  a  series  of  intermediate  time  steps  (simulating  the  manual  restarts 

required  in  KXL)  automatically  adjusting  the  length  of  each  intermediate  time  step  based 

upon  whether  or  not  the  previous  time  step  was  successful.  This  is  facilitated  by  the  ion 

current  control  and  allows  the  code  to  run  unattended  as  far  as  is  possible  without 

requiring  operator  intervention. 

*  * 

IV.  SAMPLE  CODE  RUN  FOR  EKXL  V4.1 

Figure  7  shows  a  sample  input  file  for  EKXL  version  4.1  as  it  is  described  above.  This 
calculation  tan  in  automatic  mode  with  fixed  ion  current  for  a  15  keV  well  in  a  92  cm 
device.  The  name  of  the  log  file,  a  sample  of  which  is  shown  in  Figures  8  a— b,  and  the 
dataset  file,  a  sample  of  which  is  shown  in  Figure  9,  are  taken  from  the  name  of  the  input 
file,  with  the  extension  indicating  the  nature  of  the  data.  For  example,  the  input  file  in 
Figure  7  is  named  BREM2b.INP.  The  log  file  in  Figure  8  is  named  BREM2b.LOG  and  the 
dataset  file  in  Figure  9  is  named  BREM2b.OUT.  The  last  character  of  the  file  name  is 
used  to  indicate  the  number  of  the  run.  In  this  example,  the  "b"  indicates  that  it  is  the 
second  run  and  the  "Initializ^on  file"  is  set  to  "YES"  in  the  input  file.  This  calculation 
started  where  BREM2a.INP  ended,  and  this  is  verified  in  the  beginning  of  the  log  file 
(Figure  8  a)  where  it  states  that  initialization  file  BREM2.TMP  was  read. 

The  log  file  prints  run  diagnostics  of  the  numerical  calculation^  to  help  the  user  in 
determining  if  and  when  problems  occurred  in  the  calculation.  Figure  7  a  shows  the 
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beginning  of  the  log  file,  while  Figuie  7  b  shows  the  end.  As  seen  in  the  figures,  the  log  file 
also  contains  the  transit  times  (Tau),  the  electron  confinement  time,  the  ion  upscattering 
time  (if  this  calculation  is  being  used),  the  radius  of  the  wiffle  ball,  <r^>,  the  actual  ion 
and  electron  currents  being  used  in  the  code,  the  fusion  rates  for  DT,  D  He  and  p  B,  the 
plasma  gain  for  all  four  fuels  and  the  bremmstrahlung  power  as  calculated  from  the  DD 
density.  The  end  of  the  file  indicates  that  the  snapshot  was  saved  and  lists  the  files  to 
which  output  was  written. 

Figure  9  shows  only  the  first  portion  of  an  output  dataset  file.  This  file  contains  the 
potential  and  ion  and  electron  densities  at  each  grid  point  for  every  successful  snapshot  in 
the  run.  The  output  is  formatted  and  comma-separated  which  allows  for  easy  import  into 
spreadsheet  or  graphics  software.  A  sample  plot  of  this  output  data  is  shown  in  Figure  10. 
The  beam  profiles  from  the  input  file  have  been  printed  across  the  top  of  this  chart  and  the 
fusion  rates  have  been  imported  from  the  log  file. 

Hundreds  of  cases  have  been  run  with  the  EKXL  code  modelling  many  areas  of  interest 
from  basic  experiments,  to  cunent  SCIF  parameters,  to  reactor  regimes.  Results  of  these 
code  runs  have  been  used  extensively  to  better  understand  the  physics  involved  in 
Poly  well  devices  and  also  to  improve  the  modelling,  especially  as  it  is  applied  in  the 
code.  A  compilation  of  results  of  specific  interest  to  the  SCIF  experiment  are  presented  in 
some  detail  in  another  EMC2  technical  report  on  these  computer  simulations  The 
potential  of  this  code,  howev||^  has  not  been  fully  realized  and  there  still  remain  many 
areas  of  interest  that  have  not  been  explored. 
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was  generated  in  SuperCalc  v5.0. 
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Figure  1. 

Maximum  Operational  Levels 
for  Successful  Vlasov-Poisson  "1-D"  Code  Runs 
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figure  2.  a. 
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Core  density  vs.  electron  carrent  for  a  fixed  electron  confinement.  Slope  is 
linear  for  currents  greater  1000  A. 
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Fignze  6.  Formulas  used  in  the  minor/wiffle  modd  for  Gj  calculation. 


Log  Filename:  (program  uses  XNP  filename) 

- -f 


Dataset  Name:  (program  uses  ZNP  filename)/  Initialization  file:  YES, 

- - ( -  j . 

Operating  mode:  AUTOMATIC  ,  end  tiise«1.000e-03  sec,  with  3  snapshots/dt, 


- - (  J - [  SB  ] - 

Fixed  Outer  Potential:  0.e3  Volts  at  r«0.92  m,  Ncusps*  8MR,  8WB, 

Number  of  Beam  Populations:  2,  with  repeller  effic.  :  0.9000, 

- 1 - +=_+ - j — +====+_ 

Name  Ig  (A]  EinjteV)  dE(eV]  dEperp[eV)  B(G),Gi  V0trap[eV] 

1)  Elec  18.0e03  14.5e03  0.S0e03  4.35e03  3.00e03  '20.0e03 

2)  Deut  18.0e03  lO.OOeOO  l.OOeOO  1.71e00  leJO  -11.0e03 


■+ - —  [  o== 


.J  — [SSSSSS.J 


Number  of  Thermal  Populations:  0, 

- - .fSB^f- 

Name  e[Volta]  Temp(eV] 

0)  Elec  l.el6  100. e3  20. e3 

-  — ---fSB-f-  — -f  BBSSSSS-f  — -•fBBBBBBB4'-— f  BSSSSSB^ 

Deuterium  Neutral  Density:  0.el8  /m**3,  aij:  1  ICC  mass  ratio. 

Auxiliary  Potential:  14.9e3  Volts  for  e's,  with  exponent>3.  , 

—  —  —  — ->BSBBSBS'f  —  —  —  — — -^BSBB^. 

Fractional  velocity  spread:  0,  kl  Radius  Factor:  2.0, 

- J„] - (B.—  )- 

Grid:  200  points  with  adaptive  fraction  of  0.8  and  tolerance  l.e-5  , 

Transit  Time  Correction  Factor  (Y/N):  Y,  ka  <radius>:  0.83. 

- - - 1— «]- 

Modified  B  Field  Shape:  Y,  Modified  Gj:  Y,  Adjust  Ion  Current:  N. 


Bremmstrahlung  <radius>:  1.00,  Ion  Fusion  Losses:  N 

- - - 


Flguxe  7.  Sample  input  file  foi  an  EKXL  y4.1  code  run. 


FignieS.  a 


Operating  mode  is  AUTOMATIC 
Auxiliary  potential  option  is  on. 

brein2 .  tmp 


Initialization  file  read: 


iter  ■  1>  dumax  ■  3.21B-07,  step  ■  l.OOE+00 

Electron  confinement  timet  6.389E-07 
Hiffle  radius  rb:  1.936E-01 


1) 

time 

m 

l.OOOE-06, 

determ 

■ 

l.OOOE-l-20 

Tau 

m 

6.100E-08 

3.019E-06 

2) 

time 

m 

1.263E-06, 

determ 

m 

l.OOOE-l-20 

Tau 

m 

6.100E-08 

3.027E-06 

3) 

time 

m 

1.395E-06, 

determ 

m 

l.OOOE+20 

Tau 

m 

6.100E-08 

3.019E-06 

4) 

time 

m 

1.972E-06, 

determ 

m 

l.OOOB+20 

Tau 

m 

6.100E-08 

3.022E-06 

5) 

time 

m 

2.053E-06, 

determ 

m 

l.OOCS+20 

Tau 

m 

6.100E-08 

3.021E-06 

6) 

time 

m 

1.527E-06, 

determ 

m 

1.000E420 

Tau 

s 

6.100E-08 

3.021E-06 

Figure  8.  b 


iter  ■  1,  dumax  ■  6.82E-07,  step  >  l.OOE+00 

Electron  confinement  timet  6.996E-06 
Wiffle  radius  rb:  2.692E-01 


787) 

time 

■ 

1.636E-04, 

determ 

■ 

l.OOOE-i-20 

Tau 

m 

6.113B-08 

2.911E-06 

788) 

time 

m 

1.663E-04, 

determ 

m 

1.000E4-20 

Tau 

m 

6.113E-08 

2.911E-06 

789) 

time 

m 

1.676E-04, 

determ 

m 

l.OOOE+20 

Tau 

m 

6.113E-08 

2.911E-06 

790) 

time 

m 

1.736E-04, 

determ 

m 

l.OOOE+20 

Tau 

m 

6.113E-08 

2.911E-06 

791) 

time 

m 

1.744E-04, 

determ 

m 

1.000E-»-20 

Tau 

m 

6.113E-08 

2.911E-06 

792) 

time 

m 

1.690E-04, 

determ 

m 

l.OOOE-i-20 

Tau 

m 

6.113E-08 

2.911E-06 

iter  ■ 

1>  dumax  ••  4.24E-07,  step 

m 

l.OOE+00 

Electron  confinement  timet  6.995E>06 
Wiffle  radius  rb:  2.692E-01 
t(upscatter)  ■  O.OOOE-t-00  O.OOOE+00 
Gun  current  1  is:  1.800E-f04 
Gun  current  2  is:  1.800E4'04 
Gun  current  1  is:  1.800E':-04 
Gun  current  2  is:  1.800E-:'04 

Fusion  Rate:  DT-1.313E+14  DHe3-1.951E-*-10  DHe3  n-1.643E+ll  pBll-5.434E-«-02 
Gplasmat  DO-6.624E-09  DT-1.417E-06  DHe3»9.550E-10  pBll-2.900E-18 
Pbrem:  l.OSBE-t-Ol  . 

Gicc-  9.1301-06  No/Nc  •  9.S23B+01  (<No>/Nc)*2-  1.378B-f03 

###  SNAPSHOT  SAVED  ### 

Restart  file  written:  brem2.tmp 

Print  files  written:  brem2b.out 

brem2b.tim 


Figure  8.  a.  Beranning  of  the  log  file  created  by  EKXL  for  the  run  in  Figure  7. 
b.  End  of  the  log  file  created  by  EKj^  for  the  run  in  Figure  7. 


start  date:  07/16/91  Start  time:  17:29:20.81 

Dataset :  brem2b 

Time:  1.74E-04  sec 

Number  of  Grid  Points:  200 

Deut  Beam  Neutron  Count  rate  (#/sec):  DD«1.479E'i-12,  DNBO.OOOE-t-OO 


R  (m)  Pot  (V) 

0.000000,  -1.453E-«-04, 
0.000434,  -1.453E+04, 
0.000868,  -1.453E+04, 
0.001302,  -1.453E+04, 
0.001737,  -1.4S3E+04, 
0.002171,  -1.453E+04, 
0.002605,  -1.4S3E-*-04, 
0.003039,  -1.453E+04, 
0.003473,  -1.4S3E+04, 
0.003907,  -1.4S3E<t-04, 
0.004341,  -1.453E+04, 
0.004776,  -1.453E^04, 
0.005210,  -1.453E+04, 
0.005644,  -1.453E+04, 
0.006078,  -1.453E+04, 
0.006512,  -1.453E-«-04, 
0.006946,  -1.4S3E-«'04, 
0.007380,  •1.453Et'04, 
0.007815,  -1.453E+04, 
0.008249,  -1.453E+04, 
0.008683,  -1.453E+04, 
0.009122,  -1.453E+04, 
0.009574,  -1.453E+04, 
0.010045,  -1.457E-^04, 
0.010543,  -1.464E-«-04, 
0.011077,  -1.468E+04, 
0.011655,  -1.472E-I-04, 


DD  (n/scin3)  Lambda  (m) 
1.718E+17,  O.OOOE+00, 
1.718E+17,  O.OOOE+00, 
1.718E+17,  O.OOOE+00, 
1.719E+17,  O.OOOE+00, 
1.720E+17,  O.OOOE+00, 
1.722E+17,  O.OOOE+00, 
1.724E+17,  O.OOOE+00, 
1.727E+17,  O.OOOE+00, 
1.731E+17,  O.OOOE+00, 
1.736E+17,  O.OOOE+00, 
1.744E+17,  O.OOOB+00, 
1.7S3E+17,  O.OOOE+00, 
1.765E+17,  O.OOOE+00, 
1.780E+17,  O.OOOE+00, 
1.799E+17,  O.OOOE+00, 
1.821E+17,  O.OOOE+00, 
1.845E-«-17,  O.OOOE-t-OO, 
1.870E+17,  O.OOOE+00, 
1.893B-«-17,  O.OOOE-i'OO, 
1.908E+17,  O.OOOE+00, 
1.909E+17,  O.OOOE+00, 
1.885E+17,  O.OOOE+00, 
1.820E-»-17,  O.OOOE+00, 
1.313E417,  O.OOOE+00, 
7.e05E-«-16,  O.OOOE+00, 
5.492E-t-16,  O.OOOE-t-00, 
4.027E-I-16,  O.OOOB-t-00, 


Nb  (#/m**3) 

Elec  Deut 

6.601E+20,  6.601E-t-20 

6.601E-)-20,  6.601E-»-20 

6.601E+20,  6.601E-t-20 

6.601B+20,  6.601E-I-20 

6.601E'i-20,  6.601E+20 

6.601E-«-20,  6.601E+20 

6.601B+20,  6.60lE-t-20 

6.601E+20,  6.601E+20 

6.601E-I-20,  6.601E+20 

6.60lE-^20,  6.601E-I-20 

6.601E-<-20,  6.601E-«-20 

6.601E')-20,  6.60lE-«-20 

6.601E-i'20,  6.601E+20 

6.601B-*-20,  6.601E+20 

6.601E+20,  6.601E+20 

6.601B‘t'20,  6.601E-t-20 

6.601B420,  6.601E’»-20 

6.601B4-20,  6.601B-i-20 

6.601E+20,  6.60lE-^20 

6.601E+20,  6.601E+20 

6.601E+20,  6.601E-^20 

6.601E-*-20,  6.601E+20 

6.601E'«-20,  6.601E-i-20 

5.760E+20,  5.760E+20 

4.414E<i'20,  4.414E-f20 

3.689E't'20,  3.689E-I-20 

3.150E-i-20,  3.150E+20 


Rgnzc  9.  Beginning  of  the  output  dataset  file  generated  by  the  EKXL  run  in  Figure  7. 
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Flgnze  10  Sample  plot  of  the  data  in  the  EKXL  output  file,  as  shown  in  Figure  9.  This 
graph  was  generated  in  SuperCalc  v5.0. 


